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^ . Abstract 

D , Limiting fragmentation in proton-proton, deuteron-nucleus and nucleus- 

■ nucleus collisions is analyzed in the framework of the Balitsky-Kovchegov 

equation in high energy QCD. Good agreement with experimental data 
is obtained for a wide range of energies. Further detailed tests of lim- 
' iting fragmentation at RHIC and the LHC will provide insight into the 

evolution equations for high energy QCD. 

1 Introduction 



The hypothesis of limiting fragmentation in high energy hadron-hadron colli- 
sions was suggested nearly four decades ago [1]. This hypothesis states that 
the produced particles, in the rest frame of one of the colliding hadrons, will 
approach a limiting distribution. These universal distributions describe the mo- 
mentum distributions of the fragments of the other hadron. Central to the 
original hypothesis of the limiting fragmentation in Rcf. [1] was the assumption 
that the total hadronic cross sections would become constant at large center-of- 
mass energy. If this occurred, the excitation and break-up of a hadron would be 
independent of the center-of-mass energy and distributions in the fragmentation 
region would approach a limiting curve. 
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We now know that the total hadronic cross-sections are not constant at high 
energies. Instead, to the highest energies achieved, they grow slowly with the 
center-of-mass energy y/s, with favored functional forms being either a power 
law behavior <t(s) oc s a with a sa 0.08 [2] or a er(s) oc ln 2 (s) [3]. A measurement 
of the total cross-section at the LHC [4] will further help constrain the possible 
functional forms. 

Even though the cross-sections are not constant, limiting fragmentation ap- 
pears to have a wide regime of validity. It was confirmed experimentally in 
pp,pA,nA and nucleus- nucleus collisions at high energies [5-8]. More recently, 
the BRAHMS and PHOBOS experiments at the Relativistic Heavy Ion Collider 
(RHIC) at Brookhavcn National Laboratory (BNL) have performed detailed 
studies of the pseudo-rapidity distribution of the produced charged particles 
dNch/di] for a wide range (—5.4 < rj < 5.4) of pseudo-rapidities, and for sev- 
eral center-of-mass energies (y/s NN — 19.6,62.4, 130 and 200 GeV) in nucleus- 
nucleus (Au-Au and Cu-Cu) and deuteron-nucleus (d-Au) collisions. Results 
for pseudo-rapidity distributions have also been obtained over a limited kine- 
matic range in pseudo-rapidity by the STAR experiment at RHIC [9]. These 
measurements in A-A and d-A collisions were performed for several centralities. 
In addition to the d-Au and A-A data, there are pp data at y/s = 200 GeV 
and also at y/s = 410 GeV. (In the near future, data at y/s — 500 GeV may 
become available.) These measurements have opened a new and precise win- 
dow on many scaling phenomena glimpsed at lower energies. In particular, they 
have performed detailed studies of the limiting fragmentation phenomenon. The 
pseudo-rapidity distribution ^ (where rj = rj - Y hcam is the pseudo-rapidity 
shifted by the beam rapidity ybeam = hiy/s/mp) is observed to become inde- 
pendent of the center-of-mass energy y/s in the region around rj ~ 

dN ch , r dN ch , 

_( )7> V5,6) = -^ r (,,6) I (1) 

where b is the impact parameter. 

It is worth noting that this scaling is in a strong disagreement with boost 
invariant scenarios which predict a fixed fragmentation region and a broad cen- 
tral plateau growing with energy. It would therefore be desirable to understand 
the nature of hadronic interactions that lead to limiting fragmentation and the 
deviations away from it. In a recent article, Bialas and Jezabek [47], argued that 
some qualitative features of limiting fragmentation can be explained in a two- 
step model involving multiple gluon exchange between partons of the colliding 
hadrons and the subsequent radiation of hadronic clusters by the interacting 
hadrons. In this paper, we compute limiting fragmentation within the Color 
Glass Condensate (CGC) approach [14] to high energy hadronic collisions-its 
relation to the Bialas-Jezabek model will be addressed briefly later. 

In the CGC formalism, gluon production in the limiting fragmentation region 
can be described, at leading order, in the framework of fcj_-factorization. Under 
this assumption, the inclusive cross-section for produced gluons is expressed as 
the convolution of the "unintcgratcd parton distributions" in the projectile and 
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target respectively , 

<t> A ( x i> k i± = l fe iJ-l) (p B ( x 2,k 2 ± = \k 2 ±\) 5 {2) (p ± - feii. - k 2 ±) , (2) 

times a (transverse momentum dependent) vertex squared. As we shall empha- 
size later, the name "unintegrated parton distributions" though conventional is 
somewhat imprecise because the objects 4> A B that enter in this formula differ 
from the expectation value of the gluon number operator. In eq. (2), x\ and x 2 
are the longitudinal momentum fractions of the gluons probed in the projectile 
and target, namely, 



where Yb ca m = in (v / */ m jv) 1S the beam rapidity, m N is the nucleon mass, and 
p ± = ki± + k 2 ± is the transverse momentum of the produced gluon. We should 
emphasize that the kinematics here is the 2 — > 1 cikonal kinematics, which 
provides the leading contribution to gluon production in the CGC picture. 

As we will discuss further later, /c^-factorization actually works best in the 
limiting fragmentation kinematics where x\ > 0.1 and x 2 (or vice versa). 
Clearly, fc^-factorization tells us that the cross-sections depend in general on 
both y — Ibcam and y + Ybcam- However, target limiting fragmentation implies a 
dependence on solely y — Ybeam of the spectrum of produced particles integrated 
over the transverse momentum. As we shall see in section 2, such a scaling 
emerges in a straightforward manner in the fcj_-factorization framework if a) 
the typical transverse momentum in the projectile is much smaller than the 
typical transverse momentum in the target, and b) if unitarity is preserved in 
the evolution of the target with x 2 . These two conditions are naturally fulfilled 
in the CGC picture when x 2 becomes so small that the target reaches the "black 
disk" limit. More importantly, by studying how this limit is reached, one gets 
some insight into systematic deviations away from the limiting fragmentation. 

The dynamical evolution of the unintegrated distributions with x arc de- 
scribed in the CGC by the renormalization group equations called the JIMWLK 
equations [15]. These equations, more generally, describe the x evolution of n- 
point parton correlation functions. The dynamics in the evolution at high par- 
ton densities is characterized by a "saturation momentum" Q s (x). This scale 
is the typical transverse momentum of partons in the hadron or nuclear wave- 
function 2 . Partons with p± <C Q s saturate in the wave- function with occupation 
numbers of order l/a s . The implications of saturated or "black disc" distribu- 
tions for limiting fragmentation were studied previously in the CGC approach 3 

1 Here and in the following, we call the "projectile" the nucleus A which is probed at large 
xi, and "target" the nucleus B which is probed at small values of X2- Of course, our choice 
of semantics, reminiscent of what is used in fixed target experiments, is somewhat arbitrary 
in collider experiments where the lab frame and the center of mass frame are identical. 

2 Q'j is also related to the density of color charges per unit of transverse area in the hadron 
under consideration. 

3 Limiting fragmentation of protons in the CGC was studied previously in Ref. [32] albeit 
no comparisons to data were performed in this work. 
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by Jalilian-Marian [12] and in the work of Kharzecv and Levin [30] and later 
by Kharzeev, Levin and Nardi [31]. In [12], it was a priori assumed that the 
target is very dense and appears black to the dilute partons in the projectile. 
This leads to the simple formula 

^~[ Xl tf( XufJl Z) + Xl G A ( Xl ,^)] , (4) 

for the rapidity distribution of the produced hadrons in AA collisions, where 
X\f A {x\) and xiG A (xi) are the quark and gluon distributions in the projec- 
tile nucleus. With a suitable choice of a scale fi in these distributions, this 
prescription gave a reasonable description of the fragmentation region. 

In [30], the fc_L-factorization formalism was employed, together with a "satu- 
ration" ansatz for the unintegrated gluon distribution 4 . Our approach to com- 
puting the inclusive gluon production in the saturation scenario is very similar in 
spirit to the one pioneered in [30,31]. Here however, the unintegrated gluon dis- 
tributions at small x (x < 0.01) are computed using a mean field version of the 
JIMWLK equations called the Balitsky-Kovchegov (BK) equation [16,17]. This 
approximation is strictly valid in the large N c , large mass number A and high 
energy limit. However, the BK equation may have a wider range of applicability 
beyond these asymptotic limits. Remarkably, it has been shown recently that 
the BK equation lies in the same universality class as the Fisher-Kolmogorov- 
Petrovsky-Piscounov equation in statistical mechanics, which describe a wide 
range of reaction-diffusion phenomena in nature [26]. We will discuss results 
for limiting fragmentation from solutions of the BK equation for unintegrated 
distributions in both the fixed and running coupling cases. 

For x > 0.01, the unintegrated distribution, computed using the BK equa- 
tion, is smoothly matched on to a functional form that contains key features 
expected of large x parton distributions. We will also discuss what more detailed 
comparisons to the data can teach us about these distributions. 

This article is organized as follows. In section 2, we discuss inclusive gluon 
production in high energy hadron and nuclear collisions. In the kinematical 
region of interest, k± -factorization is applicable; this allows us to relate the 
distribution of produced gluons to the unintegrated gluon distributions in the 
hadron wave- functions. Parton-hadron duality [19] is assumed to compare the 
distribution of produced gluons to those of hadrons. (The effect of fragmenta- 
tion functions can be significant at higher energies-a qualitative discussion is 
presented in section 4.) The evolution of the unintegrated gluon distributions 
(with x\ and xi respectively) in the projectile and target wave-functions is de- 
scribed by solutions of the BK equation. The fixed and running coupling forms 
of the BK equation and its solution are discussed in section 3. Results for gluon 
distributions obtained from numerical solutions of the BK equation are com- 
pared to data on limiting fragmentation in section 4. We compare our results 
to data from pp, deuteron-gold and AA collisions and discuss their dependence 

4 This ansatz is not the same ansatz as what one obtains in the McLerran-Vcnugopalan 
model [35,45]. 
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on the initial conditions, running coupling effects and on infrared regulators. 
Our results are summarized and open problems emphasized in the concluding 
section. 



2 Inclusive particle production 
in high energy collisions 

At very high energies, for y/s » Q s > pr, gluon production and fragmentation is 
the dominant mechanism for particle production. When the occupation number 
of partons is small in one of the colliding hadrons and large in the other (as is the 
case in proton-nucleus or nucleus- nucleus collisions in the fragmentation region 
of one of the nuclei), the inclusive multiplicity distribution of produced gluons 
can be expressed in the /c^-factorized form [29,34,45], 

dN s a a S iB 1 f d 2 k 



j ( t ) A( x uk±)(p B {X2, |p_L-fe_L|) • (5) 



dyd 2 Pl _ 2^C f (ttRI)(ttRI) p 2 ± J (2tt) 

Strictly speaking, the formula, as written here, is only valid at zero impact 
parameter and assumes that the nuclei have a uniform density in the transverse 
plane. Indeed, the functions (f> A B are defined for the entire nucleus. Here, S AB 
denotes the transverse area of the overlap region between the two nuclei, while 
nR A b are the total transverse area of the nuclei, and C F = (N 2 — 1)/2N C is 
the Casimir in the fundamental representation. The longitudinal momentum 
fractions X\ and x 2 were defined previously in eq. (3). 

The functions 4> A and <j) B are obtained from the dipole-nucleus cross-sections 
for nuclei A and B respectively, 



ttR 2 k 2 r 



(6) 



where Y = \n(l/x) and where the matrices U are adjoint Wilson lines evaluated 
in the classical color field created by a given partonic configuration of the nuclei 
A or B in the infinite momentum frame. For a nucleus moving in the — z 
direction, they are defined to be 



U{x ± ) = T+cxp 



-ig 



2 



J dz+^piz+^x^-T 



(7) 



Here the T a are the generators of the adjoint representation of SU(N C ) and 
T + denotes the "time ordering" along the z + axis. p a (z + ,x±) is a certain 
configuration of the density of color charges in the nucleus under consideration, 
and the expectation value ( ■ ■ ■ ) corresponds to the average over these color 
sources p a . In the McLerran-Venugopalan (MV) model [13], where no quantum 
evolution effects are included, the p's have a Gaussian distribution, with a 2- 
point correlator given by (p a (0)p b (x±)) = p 2 A 5 ab 5 ( - 2 \x 1 _ - y ± ), where p\ = 
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2 ^2 • This determines A B completely [34,45], since the 2-point correlator 

A 

is all we need to know for a Gaussian distribution. As we will discuss later, 
the small x quantum evolution of the correlator in eq. (6) is given by the BK 
equation. 

These distributions <j> A B , albeit very similar to the canonical unintcgrated 
gluon distributions in the hadrons, should not be confused with the latter [35,45]. 
However, at large k± (k± 3> Q s ), they coincide with the usual unintegrated 
gluon distribution and this determines their normalization 5 . 

The fcj_-factorized expression in eq. (5) is only valid for inclusive gluon 
production when one of the hadrons is dilute and the other is dense. In the 
CGC framework, this means that we keep only terms of orders 0{p A /k 2 ± ) and 
0(p™/fcj") in the amplitudes, if A and B are the dilute and dense hadrons re- 
spectively. This factorization is therefore applicable for proton-proton collisions, 
or nucleus- nucleus collisions in the fragmentation region of one of the projectiles. 
Clearly, it applies as well to proton/deuteron-gold collisions. fcj_-factorization 
breaks down in kinematic regions that do not satisfy this dilute-dense crite- 
rion. It particular, it is not a good approximation in the central rapidity region. 
Although there are some analytical attempts to address these violations of k±- 
factorization [36], these have been computed only numerically [37] thus far. For 
quark production, k± -factorization is broken already at leading order [46]. The 
magnitude of their breaking has been quantified recently [50] . Here we will con- 
sider only the fc^-factorized expression of eq. (5) with the understanding that 
this expression likely has significant corrections at central rapidities. These will 
be discussed further later in the paper. 

From eq. (5), it is easy to see how limiting fragmentation emerges in the 
limit where x-i <ii. In this situation, the typical transverse momentum k± in 
the projectile is much smaller than the typical transverse momentum \p ± — k± | 
in the target, because these are controlled by saturation scales evaluated respec- 
tively at xi and at x^. Therefore, at sufficiently high energies, it is legitimate 
to approximate <j) B (x2, \p± — fej_|) by cf) B (x2,p±)- When we integrate the gluon 
distribution over p ± , we thus obtain 

2^C F (nRl)(,Rl) J (2^ 
^J^^M)^^x l9 (^,^^Q 2 s (x 2 )). (8) 

This expression, in the stated approximation, is nearly independent of X2 and 
therefore depends only weakly on on y — Ybcam- To see this more clearly, note 

5 The unintegrated gluon distribution here is denned such that the proton gluon distribu- 
tion, to leading order, satisfies 

1 [® 2 

xG p (x,Q 2 ) = — J dl 2 ± cf) p (x,l ± ). 

This normalization is different from Ref. [35] - we have checked however that, when appro- 
priately normalized, our expression for the inclusive gluon distributions is identical to that of 
Ref. [35]. 



dy 
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that to obtain the second line in the above expression, we used eq. (6) and the 
fact that the Wilson line U is a unitary matrix. The details of the evolution 
are therefore not important to achieve limiting fragmentation, only that the 
evolution equation preserve unitarity. The residual dependence on X2 comes 
from the the upper limit ~ Qf (#2) of the integral in the second line. This 
ensures the applicability of the approximation that led to the expression in the 
second line above. The integral over k± gives the integrated gluon distribution 
in the projectile, evaluated at a resolution scale of the order of the saturation 
scale of the target. Therefore, the residual dependence on y + Ybeam arises only 
via the scale dependence of the gluon distribution of the projectile. This residual 
dependence on y + Ybeam is very weak at large X\ because it is the regime where 
Bjorkcn scaling is observed. 

The formula in cq. 8 was used previously in Ref. [12]. The nuclear gluon 
distribution here is determined by global fits to deeply inelastic scattering and 
Drell-Yan data. We note that the glue at large x is very poorly constrained 
at present [43]. The approach of Bialas and Jezabek [47] also amounts to us- 
ing a similar formula, although convoluted with a fragmentation function (see 
eqs. (1), (4) and (5) of [47] - in addition, both the parton distribution and the 
fragmentation function are assumed to be scale independent in this approach). 
We will discuss the effect of fragmentation functions later in our discussion. 

Though limiting fragmentation can be simply understood as a consequence 
of unitarity in the high energy limit, what may be more compelling are observed 
deviations from limiting fragmentation and how these vary with energy. These 
deviations may probe more deeply our understanding of the dynamics of both 
large x and small x modes in hadronic wavefunctions. In particular, in the small 
x case, it may provide further insight into the renormalization group equations 
that, while trivially preserving unitarity, demonstrate interesting pre-asymptotic 
behavior. These concerns provide the motivation for this detailed study with 
the Balitsky-Kovchegov renormalization group equation-to be discussed in the 
following section. 

3 Balitsky-Kovchegov equation 

We begin by briefly recapitulating key features of the Balitsky-Kovchegov (BK) 
equation and its solution [16,17]. Readers arc referred to recent review literature 
on the subject for a more detailed discussion [14,20]. The BK equation is a non- 
linear evolution equation in rapidity Y = ln(l/x) for the forward scattering 
amplitude of a quark-antiquark dipolc scattering off a target in the limit of 
very high center-of-mass energy ^fs. It was originally derived, within the dipole 
picture (which assumes the large N c limit) at small values of Bjorken x, by 
taking into account multiple rescatterings of the qq dipoles off a dense nuclear 
target. The BK equation for the amplitude is equivalent to the corresponding 
JIMWLK equation [15] of the Color Glass Condensate, in a mean field (large N c 
and large A) approximation where higher order dipole correlators are neglected. 
The parametrically suppressed N c and A contributions can, in principle, be 
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computed by solving the JIMWLK equation. In momentum space, the BK 
equation takes the form 6 



OY 



= a s (K®T)(k ±1 Y)-a s T z (k ± ,Y) , (9) 



where we denote a s = a s N c /ir. The operator K is the well known BFKL kernel 
in momentum space [21], whose action on the function T is given by 



' J2 



(K<*T\(h Y\- f d{k '^ ] k?T(k' ± ,Y)-klT(k ± ,Y) kjT(k ±1 Y) \ 

( )( ' )_ / ~W { + 7WrWj ■ 

(10) 

The function T(k±,Y) is the Bessel-Fourier transform of the dipole-target scat- 
tering amplitude T(r±,Y): 

+ OC 

f{k ± ,Y)= [ — Jo(k^)T(r ± ,Y) , (11) 
J r± 
o 

where r± is the size of the qq dipole and kj_ is its conjugate transverse momen- 
tum. The dipole amplitude T is defined in terms of the correlator of two Wilson 
lines of gauge fields in the target as 

T(r±,Y) = 1 - i-Tr {W{0)U(r ± )) y , (12) 

where we have assumed translation invariance in the transverse plane in order 
to set the quark transverse coordinate to 0. Here U is a Wilson line in the fun- 
damental representation, obtained by replacing in eq. (7) the adjoint generators 
T a by the generators t a in the fundamental representation. 

However, to evaluate the "unintegrated gluon distribution" defined in eq. (6), 
one needs to know the rapidity dependence of the correlator of two Wilson lines 
in the adjoint representation, (U(0)U^(r±)} . This can be obtained as follows. 
In the CGC framework, the weight functional for the color sources in generating 
functional can be expressed in the large N c and large A limit as a non-local 
Gaussian distribution of color sources. In this limit, one can obtain closed 
form expressions for expectation values of both the fundamental and adjoint 
correlators of Wilson lines 7 . One therefore has 

Tr (u{0)W{r ± )} = N c e - c ^ r(r± - Y) , 
Tr(U(0)UHr ± )) = Nj e -c A r{ r± ,Y) _ (13) 



6 Here we present the form of the equation for the case where the forward scattering am- 
plitude is independent of the impact parameter. A numerical study of the impact parameter 
dependence of the BK equation has been performed previously [24] ; it will be considered in 
future as an extension to this work. 

7 For a detailed discussion and relevant references, we refer the reader to Appendix A of 
Rcf. [46]. 
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where C A = N c is the Casimir in the adjoint representation. The function T is 
proportional to the variance of the non-local Gaussian weight functional in the 
generating functional and is therefore the same in both the fundamental and 
adjoint cases. As C A /C F =2 in the large N c limit, eqs. (13) and (12) give 

±-Tr(U(0)U\r ± )) Y = N c [l -T(r ± ,Y)] 2 . (14) 

Substituting the LHS side here by the RHS into cq. (6), we obtain 

2 2 2 

^,B(»,fc±) = - R * Nck± [ r ± dr ± J (k ± r ± ) [l - T A , B (r±, ln(l/x))] 2 . 



(15) 

We digress here to note that if instead we had used the correlator of two Wil- 
son lines in the fundamental representation in eq. (6), we would have obtained 
the following expression for the unintegrated gluon density 

2 2 2 

4 A , B (x,k ± ) = A f r ± dr ± J (k ± r ± ) [l - T A B (r ± ,\n(l/x))] . (16) 



The change from the correlator of Wilson lines in the fundamental representation 
to the adjoint one corresponds to the emission of the gluon from the triple 
Pomeron vertex itself. This emission is not prohibited by the Abramovsky- 
Gribov-Kancheli (AGK) cutting rules [38] for inclusive gluon production in high 
energy QCD. It was argued previously [39,40] that the numerical difference in 
the resulting rapidity distributions is rather small. In figure 1, we compare 
numerical solutions (to be discussed further shortly) for the unintegrated gluon 
distribution using the correlator of Wilson lines in the adjoint representation 
(cq. 15)) with the one using the correlator of Wilson lines in the fundamental 
representation (eq.16). The shape of the distribution is very similar, whereas 
the position of the peak is different. At Y = the position of the peaks differ 
by the ratio of the Ca/Cf — 2. At higher values of Y this difference increases 
due to the faster evolution of the correlator in the adjoint representation. 

The BK equation is the simplest evolution equation capturing both the lead- 
ing ln(l/:r) BFKL dynamics at moderately small values of x, as well as the 
recombination physics of high parton densities at very small values of a;. In the 
limit where the dipole scattering amplitude is small, T <C 1 , the non- linear term 
in eq. (9) can be ignored and the BK equation reduces to the BFKL equation. 
In this limit, the amplitude has the solution, 



2pa s Y 



exp (i + ^ r -27fe7> ' (17) 
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Figure 1: The unintegrated gluon distribution obtained from a) correlators in 
the adjoint representation (eq.15) (dashed lines) and b) correlators in the fun- 
damental representation (eq. 16) (solid lines). 



where w = 4 In 2 « 2.77, j3 = 28C(3) « 33.67 and p = \n(r 2 ± Q ) 2 . If we define 
the saturation condition as 

T(r ± ,Y) = ± at r ± = A , (18) 

one gets [27] 

Ql = Ql e c "= y where c « 4.84 . (19) 

A more careful analysis of the BK-equation close to the saturation boundary 
gives solution c « 4.88. In the strong field (saturation) limit [25] where T w 1, 
one finds 

T(r±, y) = 1 - « exp (-1 ln 2 KQ2)^ ; (20) 

with c w 4.84 and where k is an undetermined constant. In the large Y asymp- 
totic regime, leading and next-to-leading corrections to the form of the BK 
amplitude and the saturation scale have been computed [26,28]. 

The coupling constant a s in eq. (9) is fixed in the leading logarithmic ap- 
proximation in 1/x. This leads to a rather strong dependence of the saturation 
scale in the leading logarithmic approximation: Q s ~ x~ Xs with A s i=s 4.88 a s . 
For typical values of a s , of order 0.2-0.3, A s is at least a factor of three higher 
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than fits to the HERA and RHIC data would suggest [48] . It is therefore desir- 
able to include at least part of the next-to-leading corrections in ln(l/x) to the 
BK equation. Such corrections are known to reduce significantly the Pomeron 
intercept in the case of the BFKL equation. 

The BK equation has been solved numerically for both fixed and running 
coupling [22-24,42]. These solutions have the following features: 

• The dipolc amplitude is shown to unitarize, and the solution exhibits 
geometrical scaling. 

• The saturation scale has the behavior in eq. (19) or Q s ~ exp(-\/Aln TJx) 
in the fixed and running coupling cases respectively. 

• The infrared diffusion problem of the BFKL solution is cured by the non- 
linear term in the BK equation. 

• Ultraviolet diffusion is still present in the BK equation. However, it is 
attenuated by including running coupling effects [22,42,49]. 

In this work, we solved the BK equation numerically, in both fixed and run- 
ning coupling cases, to investigate limiting fragmentation in hadronic collisions. 
In the fixed coupling case, since realistic values of a s give too high a value of 
A s , we solved the BK equation with very small values of a s chosen to give X s 
values that are compatible with data from HERA, RHIC and hadron colliders. 
(We will see that these lower values are indeed favored by the data.) This is 
not unreasonable because it has been shown [28] that resummed next-to-leading 
order corrections to the BFKL equation in the presence of a unitary boundary 
(a la BK) give values of the saturation scale A s that are nearly independent of 
the energy, as in the LO case, but with a much smaller value of A s . 

The limited running coupling studies we performed did not give nearly as 
good agreement with the data than the fixed coupling studies with the lower 
values of a s . This is likely due to the fact that the energy dependence in this case 
is much too fast compared to the data. Further, since our results are sensitive to 
small transverse momenta, they will also depend strongly on how one regulates 
the running of as in the infrared. This requires a more detailed study than 
reported here. We shall therefore restrict ourselves to the fixed coupling case in 
discussing our results in the following section. 

As explained before, the "unintegrated gluon distributions" (j) A B are ob- 
tained from solutions to eq. (9). At small values of x, where the gluon density 
is high, eq. (9) captures the essential physics of saturation; we therefore use it 
to determine (f>(x, k±) for values of x < x with x = 0.01. For larger values of 
x (x > 0.01), we used the phenomenological extrapolation 8 , 

<i>(x,k±) = (~r) 4>(xo,kj_) , x > x Q . (21) 

8 A similar extrapolation was also used in Ref. [50] to study quark pair production. 
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The parameter (3 = 4 is fixed by QCD counting rules [51]. We checked that the 
results are not sensitive to the variation of this parameter in the range 4 — 5. 
The parameter A was varied between and 0.15 in the fits to data discussed 
in the following section. 

The initial condition for the BK evolution, that gives <j>(xo,k±), is given 
by the McLerran-Venugopalan (MV) model [13] with a fixed initial value of 
the saturation scale Qf(xo). For a gold nucleus, extrapolations from HERA 
and estimates from fits to RHIC data suggest that (Qf(xo)) 2 w 2GeV 2 . The 
saturation scale in the proton is taken to be Q 2 (xq) = (Qf(xo)) 2 (197/A) 1 / 3 . 
For comparison, we also considered initial conditions from the Golcc-Biernat 
and Wusthoff (GBW) model [48]. The values of Qf were varied in this study 
to obtain best fits to the data. 



4 Results for Limiting Fragmentation 

Experimental data are presented in terms of the measured distributions of pro- 
duced charged hadrons. Our expression in eq. (5) is for produced gluons. Ra- 
pidity distributions are dominated by low p± (p± < Q s ) particles; the detailed 
mechanism of how such soft gluons fragment to form hadrons is unknown. How- 
ever, in e + e~ collisions, several studies have been performed of hadronization; 
it is observed that at scales comparable to p± ~ Q s , the distribution of pro- 
duced hadrons mirrors that of the produced gluons. This hypothesis is known 
as "parton-hadron" duality [19] and we shall adopt it here. 9 We will discuss the 
effect of fragmentation functions shortly. 

Limiting fragmentation is often experimentally studied in terms of the mea- 
sured pseudo-rapidity 77 of produced particles. For massless particles, rj and y 
are the same, but they differ for massive particles. One obtains 



y(r],p±,m) = ^ln 



\Jm 2 +p 2 L ch 2 ?7 + p± sh?7 
\Jm 2 +p 2 L ch 2 ?7 — p± sh?7 



(22) 



We choose m to be of order of m ~ 150 MeV — 300 McV. Our expression 
for dN/dy of gluons has a logarithmic divergence at low p±; in addition to 
parton-hadron duality, for self consistency, we should regulate the corresponding 
expression for charged hadrons with the same mass as the one used in the 
conversion from y to r\. We will later discuss the sensitivity of our results to m. 

We now have all the ingredients necessary to calculate the multiplicity dis- 
tributions in the Color Glass Condensate framework. As discussed previously, 



9 This assumption was also made in rcf. [30] and is implicit in several other works. Note that 
such an assumption may be invalidated by final state interactions such at medium induced 
parton splitting, namely, "jet quenching". As we shall discuss later in detail, we observe that 
the p± spectral shapes of softer partons, albeit one would imagine them to interact more 
strongly, are closer to that of observed hadrons-than those with > Qs- This is fortuitous 
for the discussion of limiting fragmentation here since dN/ dy is dominated by soft momenta. 
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experimental data on pseudo-rapidity distributions in the projectile fragmenta- 
tion region scale with 77 — Ybeam for different energies. In eq. (5), 4> A {x\) depends 
on the difference y — Ybeam and B (x 2 ) on the sum y + Ybeam- To provide a 
sense of the x's involved, consider gold-gold collisions at RHIC. For y = Yb oa m, 
which lies in the fragmentation region, and p± <~ <~ 1 GeV , Xi s=s 1 and 
x 2 ~ 2.5 • 10~ 5 . One is therefore probing very small values of x 2 in nuclei in 
this region 10 , where the saturation scale Q s (x 2 ) is rather large. As we discussed 
previously, in this kinematic region, because of unitarity, the gluon distribution 
depends only weakly on x 2 ~ cxp(— (y + Yb cam )). This qualitatively explains 
the scaling in the limiting fragmentation region. 

We now turn to our results. To reiterate, they are obtained by a) solving 
the BK equation to obtain cq. (11), and thereby eq. (15) for the unintegrated 
distribution for x < 10~ 2 , b) using eq. (21) to determine the large x (x > 10~ 2 ) 
behavior, and c) substituting these in eq. (5) to determine rapidity distribution 
of gluons. The pseudo-rapidity distributions are determined by the transforma- 
tion in eq. (22). 

In figure 2 we show pseudo-rapidity distributions of the charged particles 
produced in nuclcon-nuclcon collisions for center of mass energies ranging from 
53 GeV to 900 GeV. We have performed calculations for two different input 
distributions at starting value of x — 0.01, namely, the GBW and MV models. 
The normalization has been treated as a free parameter. Extrapolations to large 
values of x are performed using the formula in eq. (21). Plots on the left hand 
side of figure 2 are obtained for Ao = 0.15 whereas the right hand side plots are 
done for Ao = 0.0. The different values of A s are obtained by varying a s when 
solving the BK equation. 

Note that while A s controls the growth of Q s with energy, the amplitude in 
cq. (17) has the growth rate 

2.77 

Abk = A s w 0.57 X s , 

which is significantly lower. So A s = 0.46, which gives reasonable fits (more 
on this in the next paragraph) to the pp data for the MV initial conditions, 
corresponds to Abk = 0.28, which is close to the value for the energy dependence 
of the amplitude in ncxt-to- leading order resummed BK computations [28] and 
in empirical dipole model comparisons to the HERA data [48] . 

We find that our computations are extremely sensitive to the extrapolation 
prescription to large x. This is not a surprise since we are probing the wave- 
function of the projectile at fairly large values of X\. From our analysis, we 
see that the data naively favor a non-zero value for Ao in eq. (21). The zero 
value of Ao results in the distributions which, in both the MV and GBW cases 
give reasonable fits (albeit with different normalizations) at lower energies but 
systematically become harder relative to the data as the energy is increased. To 
fit the data in the MV model up to the highest UA5 energies, a lower value of 

10 Note that x\ can be larger than unity in a nucleus and could in principle take values up 
to A, the number of nucleons. 
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Figure 2: Pseudorapidity distributions dN/d/q for charged particles from 
nucleon-nucleon collisions at UA5 energies [5] ^fs = 53, 200, 546, 900 GeV and 
PHOBOS data [8] at «Js = 200 GeV. Upper plots: initial distribution from the 
MV model, lower plots: initial distribution from the GBW model. Left panels: 
A = 0.15, right panels Aq = 0.0. 



A s than that in the GBW model is required. This is connected with the fact 
that MV model has tails which extend to larger values in k± than in the GBW 
model. As the energy is increased, the typical k± ~ Q s {x2) does as well. We 
will return to this point shortly. 

In figure 3 the same distributions are shown as a function of the r{ = 
V ~ ^bcam- The calculations for Ao = 0.15, are consistent with scaling in the 
limiting fragmentation region. There is a slight discrepancy between the cal- 
culations and the data in the mid-rapidity region. This discrepancy may be a 
hint that one is seeing violations of fcj_ formalism in that regime because k± 
formalism becomes less reliable the further one is from the dilute-dense kine- 
matics of the fragmentation regions [37,36]. This discrepancy should grow with 
increasing energy. However, our parameters are not sufficiently constrained that 
a conclusive statement can be made. 

In changing from rapidity to pseudo-rapidity distributions, one has to use a 
Jacobian with a mass parameter. We have chosen this parameter, for consis- 
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Figure 3: The same as figure 2 but plotted versus rf = rj — Ybcam to illustrate 
the region of limiting fragmentation. 



tency, to be the same as our p± cut-off for the integration over the transverse 
momentum in the formula for the gluon production. We have checked the 
sensitivity of our calculations to variations in the p± cut-off m in the range 
150 — 300 MeV. In figure 4 we show the results for the pp collisions with two 
different values of pj_. m in. = 150,300GeV. The 'dip' in the pseudo-rapidity dis- 
tribution becomes more or less pronounced when the parameter m is decreased 
or increased respectively. We note that in order to obtain a reasonable descrip- 
tion of the data we also had to adjust the other parameters (normalization and 
A.). 

In figure 5 we show the extrapolation to higher energies, in particular the 
LHC range of energies for the calculation with the GBW input. We observed 
previously that the MV initial distribution, when evolved to these higher ener- 
gies, gives a rapidity distribution which is very flat in the —5 < rj < 5 regime. 
We noted that this is because the average transverse momentum grows with the 
energy giving a significant contribution from the high k± tail of the distribution 
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Figure 4: Calculations for different values of the infrared pr cutoff compared to 
the pp data in figure 2. 



in the MV input in eq. (21). 

The effect of fragmentation functions on softening the spectra in the limiting 
fragmentation region can be simply understood by the following qualitative 
argument. The inclusive hadron distribution can be expressed as 



dN h 



d 2 p±dy 



dz dN„ 



z d 2 q±dy 



D 



g—>h 



PA 



(23) 



where -D s ^/j(z, /i 2 ) is the fragmentation function denoting the probability, at 
the scale fj?, that a gluon fragment into a hadron carrying a fraction z of its 
transverse momentum. For simplicity, we only consider here the probability 
for gluons fragmenting into the hadron. The lower limit of the integral can be 
determined from the kinematic requirement that x±_2 < 1-we obtain, 



,y-Y bB 



m N 



(24) 



16 



if z m in. were zero, the effect of including fragmentation effects would simply be 
to multiply cq. (23) by an overall constant factor. At lower energies, the typical 
value of q± is small for a fixed y— ^beam! the value of z m ; n . is quite low. However, 
as the center of mass energy is increased, the typical q± value grows slowly with 
the energy. This has the effect of raising z mm . for a fixed y — Ibeam, thereby 
lowering the value of the multiplicity in eq.(23) for that y — lb oa m- Note further 
that eq. (24) suggests that there is a kinematic bound on q± as a function of 
y — lbcam _ only very soft gluons can contribute to the inclusive multiplicity. 

In figure 6 we illustrate p± distributions obtained from the MV input com- 
pare to the UA1 data [52]. We compare the calculation with and without the 
fragmentation function. The fragmentation function has been taken from [53]. 
Clearly the "bare" MV model does not describe the data at large k± because 
it does not include fragmentation function effects which, as discussed, make the 
spectrum less flat. In contrast, because the k± spectrum of the GBW model 
dies exponentially at large k±, this "unphysical" k± behavior mimics the effect 
of fragmentation functions-sec figure 6. Hence extrapolations of this model, as 
shown in figure 5 give a more reasonable looking result. Similar conclusions 
were reached previously in Ref. [44]. 

We next compute the pseudo-rapidity distribution in dcutcron-gold colli- 
sions. In figure 7 we show the result for the calculation compared with the dA 
data [8]. The unintegrated gluons were extracted from the pp and AA data. 
The overall shape of the distribution matches well on the deuteron side with the 
minimum-bias data. The disagreement on the nuclear fragmentation side is easy 
to understand since, as mentioned earlier, it requires a better implementation of 
nuclear geometry effects. Similar conclusions were reached in Ref. [31] in their 
comparisons to the RHIC deuteron-gold data. 

We now apply our considerations to central nucleus-nucleus collisions. In 
figure 8 we present fits to data on the pseudo-rapidity distributions in gold- 
gold collisions from the PHOBOS, BRAHMS and STAR collaborations. The 
data [8] are for y/1 = 19.6, 130, 200 GeV and the BRAHMS data [7] are for 
y/s — 130, 200 GeV. A reasonable description of limiting fragmentation is 
achieved in this case as well. One again has discrepancies in the central ra- 
pidity region as in the pp case. As mentioned previously, a natural explanation 
for this discrepancy is the violation of /c^-factorization. It is expected that 
these violations decrease the multiplicity in this regime [37,36] relative to ex- 
trapolations using fcj_-factorization. We find that values of Qf w 1.3 GeV for 
the saturation scale give the best fits. This value is consistent with other es- 
timates [37,30,33]. Apparently the gold-gold data are better described by the 
calculations which have Ao = 0.0. This might be related to the difference in 
the large x distributions in the proton and nucleus. Also, slightly higher values 
of X s are preferred. This variation of parameters from A A to pp case might 
be also connected with the fact that in our approach the impact parameter is 
integrated out, so that there is no detailed information on the nuclear geometry 
A more detailed calculation with the impact parameter dependence taken into 
account is left for future investigations. 

In fig. 9 we show the extrapolation of two calculations to higher energy 
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Figure 5: Predictions for higher center-of-mass energies: y/s = 2, 8,14TeV for 
GBW input model. Parameter in the large x extrapolation was set to Aq = 0.15. 



yfs = 5500 GeV. We note that the calculation within the MV model gives 
results which would violate the scaling in the limiting fragmentation region by 
approximately 25% at larger y — Ybcam- This violation is partly because of 
the effect of fragmentation functions discussed previously and in part due to 
the fact that the integrated parton distributions from the MV model do not 
obey Bjorken scaling at large values of x. In the latter case, the violations are 
proportional to \n(Q 2 s (x2)) as discussed previously. The effects of the former are 
simulated by the GBW model-the extrapolation of which, to higher energies, is 
shown by the dashed line. The band separating the two therefore suggests the 
systematic uncertainity in such an extrapolation coming from a) the choice of 
initial conditions and b) the effects of fragmentation functions which are also 
uncertain at lower transverse momenta. 

Finally, the transverse momentum distribution in the gold-gold collisions 
is presented in fig. 10. The data at ■s/a = 200GeV are from the PHOBOS 
collaboration [11]. Again, in this case, the calculation without fragmentation 
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function tends to overshoot the data significantly at high values of pt- Including 
fragmentation function effects [53] results in a better agreement with the shape 
of the p t distribution as expected. A more detailed study to extract parameters 
that give the best chi-squared values will be performed at a later date. 

5 Conclusions 

We studied here the phenomenon of limiting fragmentation in the Color Glass 
Condensate framework. In the dilute-dense (projectile-target) kinematics of 
the fragmentation regions, one can derive (in this framework) an expression 
for inclusive gluon distributions which is k± factorizable into the product of 
"unintegrated" gluon distributions in the projectile and target. From the general 
formula for gluon production(eq. (5)), limiting fragmentation is a consequence 
of two factors: 

• Unitarity of the U matrices which appear in the definition of the uninte- 
grated gluon distribution in eq. (6). 

• Bjorkcn scaling at large x\, namely, the fact that the integrated gluon 
distribution at large x, depends only on x\ and not on the scale Q s (x 2 ). 
(The residual scale dependence consequently leads to the dependence on 
the total center-of-mass energy.) 

Deviations from the limiting curve at experimentally accessible energies are 
very interesting because they can potentially teach us about how parton distri- 
butions evolve at high energies. In the CGC framework, the Balitsky-Kovchegov 
equation determines the evolution of the unintegrated parton distributions with 
energy from an initial scale in x chosen to be x§ = 0.01. This choice of scale is 
inspired by model comparisons to the HERA data. 

We compared our results to data on limiting fragmentation from pp collisions 
at various experimental facilities over a wide range of collider energies, and to 
collider data from RHIC for deuteron-gold and gold-gold collisions. We obtained 
results for two different models of initial conditions at x > xq] the McLerran- 
Venugopalan model (MV) and the Golec-Biernat-Wusthoff (GBW) model. In 
addition to the two parameters in the initial conditions (A s and Ao), we also 
studied the sensitivity of our results to an infrared momentum cut-off m (chosen 
to be the same value in the y —* i) conversion for hadrons). 

We found reasonable agreement for this wide range of collider data for the 
limited set of parameters. Clearly these can be fine tuned by introducing fur- 
ther details about nuclear geometry. That would introduce further parameters 
but on the other hand there is more data for different centrality cuts as well-we 
leave these detailed comparisons for future studies. In addition, an important 
effect, which improves agreement with data, is to account for the fragmentation 
of gluons in hadrons. In particular, the MV model, which has the right lead- 
ing order large k± behavior, but no fragmentation effects, is much harder than 
the data. The latter falls as a much higher power of k±. Since even rapidity 
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distributions at higher energies data are more sensitive to larger kj_, we expect 
this discrepancy to show up in our studies of limiting fragmentation and indeed 
it does. We noticed that taking this into account lead to much more plausible 
extrapolations of fits of existing data to LHC energies. This "gluon fragmen- 
tation" contribution also suggests that the "flat" deviations that we found (for 
fits with MV initial conditions) for X s = 0.46 (Abk = 0.28) are diminished. 

Clearly, the procedure employed in this paper needs further improvements. 
One of them is the impact parameter dependence of the unintegrated gluon 
distribution functions. We discussed briefly "gluon fragmentation" effects (with 
different fragmentation function sets) which need to be taken into account. Fur- 
thermore, large x distributions also need to be better constrained and consis- 
tency with computations of other final states in the CGC framework established. 
Finally, the factorization formula used in this paper involves only gluons; for 
more realistic estimates, one should also include quark distributions in this 
framework. 
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Figure 6: p_L distribution from cq. (5) with MV (full squares) and GBW (full 
triangles) initial conditions. The MV initial condition- with the fragmentation 
function included-is denoted by the open squares. The distribution is aver- 
aged over the rapidity region y = 0.0 — 2.5, to compare with data (in 200 
GeV/nucleon proton-antiproton collisions in the same pseudo-rapidity range) 
on charged hadron p± distributions from the UA1 collaboration: full circles. 
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d-Au Vs = 200GeV 




Figure 7: Comparison of computations to minimum bias deuteron-gold data at 
RHIC [8]. Calculation done using the MV model as a starting distribution for 
the evolution. 
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Figure 8: Pseudorapidity distributions normalized by the number of participants 
for charged hadrons in gold-gold collisions from the PHOBOS collaboration at 
energies 19.6, 130, 200 GeV (filled triangles , squares and circles), BRAHMS col- 
laboration at energies 130, 200 GeV (open squares and circles) . The data from 
the STAR collaboration at energy 62.4 GeV (open triangles) are not visible on 
this plot but can be seen more clearly in fig. 9. Upper plots: initial distributions 
from the MV model; lower plots: initial distributions from the GBW model. 
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Figure 9: Extrapolation of calculations for gold-gold collisions shown in fig. 8 
to the LHC energy y/s = 5500 GeV/ nucleon. For comparison, the same data 
at lower energies are shown. (See fig. 8.) The band is an estimate of the 
systematic uncertainty of our approach. Its lower border corresponds to the 
GBW input with Ao = 0, X s — 0.46, and its upper border to the MV input with 
A = 0, A s = 0.46. 
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Figure 10: p± distribution from eq. (5) with MV (full triangles) initial con- 
ditions. The MV initial condition with the fragmentation function included 
is denoted by the full squares. The distribution is averaged over the pseudo- 
rapidity region y = 0.2—1.4, to compare with data in 200 GeV/nucleon gold-gold 
collisions (open circles). 
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